!###############################################
 module domain
! Domain and grids
!###############################################

!*******************
  contains
!*******************
 
!===================================================================
 subroutine gridgen
! Eulerian grid and (if requested) initial location of the markers
!===================================================================
   use share
   implicit none
   integer :: i,j,k,i1,i2,j1,j2,ipor,nind,ntot,nold,Nxm,l,m,q,nblock
   integer, dimension(10) :: typor,npor
   real (kind=8) :: x,z,length
   real (kind=8), dimension(10) :: lenpor,delta0,delta1
   real (kind=8), dimension(:), allocatable :: xpom,zpom,x1,z1,x2,z2
   character (len=1) :: ch
   character (len=100) :: readval
!===================================================================


!------------------------------------------------------------------------------
!                             Eulerian grid
!------------------------------------------------------------------------------
!  dimx => real. Width of the domain 
!  dimz => real. Height of the domain 
!  nx => integer. Number of internal nodes along x 
!  nz => integer. Number of internal nodes along z
!  xgrid => real, array (0:Nx+1). x-coordinates of the grid nodes (black squares)
!  zgrid => real, array (0:Nz+1). z-coordinates of the grid nodes (black squares)
!  delx => real, array (Nx+1). Grid spacings along x
!  delz => real, array (Nz+1). Grid spacings along z
!  rgx => integer. If 1, uniform x-grid. If 0 non-uniform x-grid
!  rgz => integer. If 1, uniform z-grid. If 0 non-uniform z-grid
!------------------------------------------------------------------------------

   call readinput("uniform_x_grid",readval)
   read(readval,*) rgx
   
   call readinput("uniform_z_grid",readval)
   read(readval,*) rgz
     
   dimz=1d0
   call readinput("aspect_ratio",readval)
   read(readval,*) dimx

!-------------------
!  Setup x-grid
!-------------------

!-----------------
! Uniform grid
!-----------------
   if(rgx.eq.1) then 

     call readinput("x_points",readval)
     read(readval,*) nx

     allocate(xgrid(0:nx+1),delx(nx+1))
     dx=dimx/(nx+1)
     delx=dx
     x=0d0
     do i=0,nx+1
       if(i.gt.0) x=x+dx
       xgrid(i)=x     
     end do

   else 
!------------------------------------------------------------------------------
! Non-uniform grid
!  ipor => integer. Number of portions of dimx or dimz with different spacings
!  typor => integer, array(10). typor(i) is the type of discretization chosen 
!           for the i-th portion:
!           typor(i)=1 -> uniform discretization
!           typor(i)=2 -> non-uniform discretization with prescribed size of
!                         the first and last interval
!  npor => integer, array(10). npor(i) is the number of intervals that are  
!                              used to cover the i-th portion
!  lenpor => real, array(10). lenpor(i) is the length of the i-th portion
!  delta0 => real, array(10). delta0(i) is the size of the first interval 
!                             of the i-th portion 
!                             (to be specified only if typor(i)=2)
!  delta1 => real, array(10). delta1(i) is the size of the last interval 
!                             of the i-th portion 
!                             (to be specifie9 only if typor(i)=2)
!------------------------------------------------------------------------------
     call readinput("x_sections",readval)
     read(readval,*) ipor

     ntot=0
     nold=0
     do i=1,ipor
       write (ch, '(I1)') i
       call readinput("x_type_"//ch,readval)
       read(readval,*) typor(i)
       call readinput("x_len_"//ch,readval)
       read(readval,*) lenpor(i)
       call readinput("x_points_"//ch,readval)
       read(readval,*) npor(i)
       if(typor(i).eq.2) then
         call readinput("x_delta1_"//ch,readval)
         read(readval,*) delta0(i)
         call readinput("x_delta2_"//ch,readval)
         read(readval,*) delta1(i)
       endif
     end do

     nind=0
     length=0d0
     do i=1,ipor
       nind=nind+npor(i)
       length=length+lenpor(i)
     end do
     if(length.ne.dimx) then
       print *, 'length=',length,'dimx=',dimx
       print *, 'The sum of domain portions along x does not cover the whole dimx'
     endif

     nx=nind-1
     allocate(xpom(0:nind),xgrid(0:nx+1),delx(nx+1))
     xpom=0d0; xgrid=0d0; delx=0d0

     do i=1,ipor
       nold=ntot 
       xpom=0d0
       call algrid(ntot,typor(i),npor(i),lenpor(i),delta0(i),delta1(i),xpom)
       if(typor(i).eq.2) then
         xpom=xpom*lenpor(i)
       endif     
       if(i.gt.1) then
         do j=2,i
           xpom=xpom+lenpor(j-1)
         end do
       endif
       do j=1,npor(i)       
         xgrid(nold+j)=xpom(j)
       end do
     end do

     do i=1,nx+1
       delx(i)=xgrid(i)-xgrid(i-1)
     end do

   endif

!-------------------
!  Setup z-grid
!-------------------

!-----------------
! Uniform grid
!-----------------
   if(rgz.eq.1) then 

     call readinput("z_points",readval)
     read(readval,*) nz
     allocate(zgrid(0:nz+1),delz(nz+1))
     dz=dimz/(nz+1)
     delz=dz
     z=0d0
     do j=0,nz+1
       if(j.gt.0) z=z+dz
       zgrid(j)=z
     end do

   else 
!--------------------
! Non-uniform grid
!--------------------
     call readinput("z_sections",readval)
     read(readval,*) ipor

     ntot=0
     nold=0
     do i=1,ipor
       write (ch, '(I1)') i
       call readinput("z_type_"//ch,readval)
       read(readval,*) typor(i)
       call readinput("z_len_"//ch,readval)
       read(readval,*) lenpor(i)
       call readinput("z_points_"//ch,readval)
       read(readval,*) npor(i)
       if(typor(i).eq.2) then
         call readinput("z_delta1_"//ch,readval)
         read(readval,*) delta0(i)
         call readinput("z_delta2_"//ch,readval)
         read(readval,*) delta1(i)
       endif
     end do
   
     nind=0
     length=0d0
     do i=1,ipor
       nind=nind+npor(i)
       length=length+lenpor(i)
     end do
     if(length.ne.dimz) then
       print *, 'length=',length,'dimz=',dimz
       print *, 'The sum of domain portions along z does not cover the whole dimz'
     endif

     nz=nind-1
     allocate(zpom(0:nind),zgrid(0:nz+1),delz(nz+1))
     zgrid=0d0; delz=0d0

     do i=1,ipor
       nold=ntot 
       zpom=0d0
       call algrid(ntot,typor(i),npor(i),lenpor(i),delta0(i),delta1(i),zpom)
       if(typor(i).eq.2) then
         zpom=zpom*lenpor(i)
       endif     
       if(i.gt.1) then
         do j=2,i
           zpom=zpom+lenpor(j-1)
         end do
       endif
       do j=1,npor(i)       
         zgrid(nold+j)=zpom(j)
       end do
     end do

     do i=1,nz+1
       delz(i)=zgrid(i)-zgrid(i-1)
     end do

   endif

!---------------------------------------------------------------------
! Store the grid on file in a format suitable for plotting in GMT
!---------------------------------------------------------------------
   open(1,file='mesh.dat')
   do i=0,nx+1
     write(1,'(E14.8,2X,E14.8)') xgrid(i), zgrid(0)
     write(1,'(E14.8,2X,E14.8)') xgrid(i), zgrid(nz+1)
     write(1,'(A1)') '>'
   end do
   do j=0,nz+1
     write(1,'(E14.8,2X,E14.8)') xgrid(0), zgrid(j)
     write(1,'(E14.8,2X,E14.8)') xgrid(nx+1), zgrid(j)
     write(1,'(A1)') '>'
   end do

!---------------
!  Markers
!----------------------------------------------------------------------------
!  Nxm => integer. Number of markers in each grid cell
!  Nmk => integer. Total number of markers
!  xmk => real, array (Nmk). x-coordinates of each marker
!  zmk => real, array (Nmk). z-coordinates of each marker
!----------------------------------------------------------------------------
   call readinput("markers",readval)
   read(readval,*) markers

! Markers randomly distributed within each grid cell and all over the domain
   if(markers.eq.1) then

     call readinput("cell_markers",readval)
     read(readval,*) Nxm

     Nmk=Nxm*(nx+1)*(nz+1)
     allocate(xmk(Nmk),zmk(Nmk),mktype(Nmk))
     xmk=0d0; zmk=0d0

     l=1
     m=1
     z=0d0
     do j=0,nz
       z=zgrid(j)
   
       x=0d0
       do i=0,nx
         x=xgrid(i)

         do q=1,Nxm
           m=l+3
           xmk(l)=x+ran(l)*(xgrid(i+1)-xgrid(i))
           zmk(l)=z+ran(m)*(zgrid(j+1)-zgrid(j))
           l=l+1
         end do

       end do
     end do

     if((debug.eq.1)) then 
       print *, 'Number of markers', Nmk
       print *, 'Number of markers in each cell', Nxm
     endif

! Single markers distributed at cell centers contained in specified rectangular blocks
   elseif(markers.eq.2) then
   
     Nmk=(nx+1)*(nz+1)
     allocate(xmk(Nmk),zmk(Nmk))
     xmk=0d0
     zmk=0d0

     call readinput("markers_blocks",readval)
     read(readval,*) nblock !  number of rectangular blocks where to put tracing particles

     allocate(x1(nblock),z1(nblock),x2(nblock),z2(nblock))

! Read coordinates of bottom-left and top-right angles of each rectangle that will contain markers
     do k=1,nblock
       write (ch, '(I1)') k

       call readinput("xblock_BL_"//ch,readval)
       read(readval,*) x1(k)
       call readinput("zblock_BL_"//ch,readval)
       read(readval,*) z1(k)

       call readinput("xblock_TR_"//ch,readval)
       read(readval,*) x2(k)
       call readinput("zblock_TR_"//ch,readval)
       read(readval,*) z2(k)

     end do

     l=0
     do k=1,nblock

       i1=bisec(xgrid,x1(k))
       i2=bisec(xgrid,x2(k))

       j1=bisec(zgrid,z1(k))
       j2=bisec(zgrid,z2(k))
       
       do i=i1,i2
         do j=j1,j2
           l=l+1
           xmk(l)=xgrid(i) + delx(i+1)/2d0
           zmk(l)=zgrid(j) + delz(j+1)/2d0
         end do
       end do

     end do

   endif

 end subroutine gridgen 


!========================================================================================
 subroutine algrid(ntot,typ,n,length,d0,d1,grid)
! Algebraic grid generation 
! for details, see Ch. 8 of Thompson et al. 'Numerical grid generation'. Elsevier, 1985
!========================================================================================
   implicit none
   integer :: n,i,typ
   integer, intent(inout) :: ntot
   real (kind=8) :: length,delta,deltaold,alfa,beta,s01,x,u,     &
                    d0,d1,dx0,dx1,dx,ds0,ds1,treq,dtreq,diff
   integer, parameter :: nmax=100000
   real (kind=8), parameter :: eps=1d-6
   real (kind=8), dimension(0:) :: grid
!=================================================

! Trascendental equation to determine delta and its derivative
   treq(delta)=sinh(delta)/delta-beta
   dtreq(delta)=(delta*cosh(delta)-sinh(delta))/delta**2

! Stretching function:   
   u(x)=0.5d0*(1d0+tanh(delta*(x-0.5d0))/tanh(delta/2))
   s01(x)=u(x)/(alfa+(1d0-alfa)*u(x))
!===========================================================   
 

!--------------------------
! Uniform grid
!--------------------------
   if(typ.eq.1) then

     dx=length/n
     do i=1,n
       grid(i)=grid(i-1)+dx
     end do
     ntot=ntot+n

!--------------------------
! Algebraic grid
!--------------------------
   elseif(typ.eq.2) then

     dx=1d0/n
     dx0=d0/length
     dx1=d1/length

     ds0=dx0/dx
     ds1=dx1/dx
     alfa=dsqrt(ds1/ds0)
     beta=1d0/dsqrt(ds1*ds0)

! Determine the proper delta iteratively (Newton method)  
     delta=1d0
     do i=1,nmax
       deltaold=delta
       delta=deltaold-treq(delta)/dtreq(delta) 
       diff=dabs(delta-deltaold)
       if(diff.le.eps) exit
     end do

! Find the grid points 
     x=0d0
     do i=1,n
       x=x+dx
       grid(i)=s01(x)
     end do
     ntot=ntot+n

   endif  

 end subroutine algrid



 end module domain
!###############################################



